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ABSTRACT 



n We investigate the couphng between rock-size sohds and gas during the for- 

P^ mation of gas giant planets by disk fragmentation in the outer regions of massive 

^ disks. In this study, we use three-dimensional radiative hydrodynamics simula- 

I tions and model solids as a spatial distribution of particles. We assume that half 

^ of the total solid fraction is in small grains and half in large solids. The former 

^ are perfectly entrained with the gas and set the opacity in the disk, while the 

I— ' latter are allowed to respond to gas drag forces, with the back reaction on the 

psj gas taken into account. To explore the maximum effects of gas-solid interactions, 

^ we first consider lOcm-size particles. We then compare these results to a simula- 

CN tion with 1 km-size particles, which explores the low-drag regime. We show that 

P^ (1) disk instability planets have the potential to form large cores due to aerody- 

l^ namic capturing of rock-size solids in spiral arms before fragmentation; (2) that 

O temporary clumps can concentrate tens of M^ of solids in very localized regions 

,-H before clump disruption; (3) that the formation of permanent clumps, even in 

> the outer disk, is dependent on the grain-size distribution, i.e., the opacity; (4) 

K> that nonaxisymmetric structure in the disk can create disk regions that have a 

solids-to-gas ratio greater than unity; (5) that the solid distribution may affect 

the fragmentation process; (6) that proto-gas giants and proto-brown dwarfs can 

start as differentiated objects prior to the H2 collapse phase; (7) that spiral arms 

in a gravitationally unstable disk are able to stop the inward drift of rock-size 

solids, even redistributing them to larger radii; and, (8) that large solids can 

form spiral arms that are offset from the gaseous spiral arms. We conclude that 

planet embryo formation can be strongly affected by the growth of solids during 

the earliest stages of disk accretion. 



1. Introduction 

Core accretion (e.g., Pollack et al. 1996; Hubickyj et al. 2005) and direct formation by 
disk instability (e.g, Cameron 1978; Boss 1997, 1998) are both viable modes of gas giant 
planet formation, and both may be required to explain the diversity of observed planetary 
systems. For disk instability to lead to fragmentation, the Toomre (1964) Q parameter must 
be rapidly lowered in the outer disk, which is where the disk is most likely to have both a 
low Q and short cooling times relative to the local orbital time (e.g., Rafikov 2007, 2009; 
Nero & Bjorkman 2009; see Durisen et al. 2007 for a review). Although disk instability 
is quite sensitive to the details of radiative transfer approximations and to opacity (e.g., 
Johnson & Gammie 2003; Cai et al. 2006, 2008, 2009; Boley et al. 2007b; Mayer et al. 2007; 
Cossins et al. 2009), hydrodynamics simulations of extended, outer disks demonstrate that 
fragmentation is possible even under nonisothermal conditions (e.g., Voroybov & Basu 2006, 
2010; Stamatellos et al. 2007; Stamatellos & Whitworth. 2009; Stamatellos et al. 2009; Boley 
2009; Forgan et al. 2009; Hayfield et al. 2010). If fragmentation is to occur, it is likely to 
happen during the earliest stages of disk evolution, when mass infall onto the disk can 
proceed faster than slow, diffusion-like mass transport (e.g., Boley 2009; hereafter B2009). 
This mass loading will likely lead to the development of global spiral modes, nonlocal mass 
transport, and in some cases, fragmentation. This led B2009 to suggest two modes of gas 
giant planet formation, with disk instability representing the early mode, operating in newly- 
formed systems, and with core accretion representing the dominant mode during the evolved 
stages of disk evolution. An example of a system with disk instability planets may very well 
be HR8799 (Marois et al. 2008), as other formation scenarios seem to be unlikely (Dodson- 
Robinson et al. 2009). 

Several observational signatures have been proposed to test whether disk instability 
plays a role in planet formation. B2009 suggested that a bimodal population of gas giants 
would reflect that both formation channels are possible and that the ratio of wide orbit gas 
giants to the total number of planets in a system should increase with decreasing metallicity. 
The latter expectation reflects that core accretion should be more sensitive to changes in 
metallicity than disk instability. If, instead, wide orbit gas giants are produced mainly by 
scattering, Dodson-Robinson et al. (2009) pointed out that the frequency of systems with 
gas giants on wide orbits should be dependent on system age. As more systems with wide 
orbit gas giants are discovered, these predictions will be tested against observations. 

In addition to population trends, metallicity and core size may be a way to distinguish 
between formation mechanisms. For example. Helled & Bodenheimer (2010) use a planetary 
evolution code to calculate contraction times for isolated clumps, and follow planetesimal 
capture of large bodies (> 1 km) until molecular hydrogen dissociates in the clump's core. 



This dissociation causes the clump to collapse rapidly to about a Juptier radius, which 
significantly reduces the cross section of the young planet and marginalizes planetesimal 
capture. They find that planetesimal capture can only contribute negligible to modest en- 
hancement of solids for planets at large radii, concluding that planets that form on wide 
orbits by disk instability should have nearly the same composition as their host star. If this 
is correct, metallicity could be a valuable tool in determining formation modes for wide orbit 
gas giants, as gas giant planets that form by core accretion are expected to be metal rich 
(e.g.. Pollack et al. 1996). However, a limitation of the Helled & Bodenheimer calculation is 
that their simulations are not within the context of a gravitationally unstable disk. When 
clumps form, they tend to form near corotation or at intersections of spiral arms (Durisen 
et al. 2008), yielding initial fragment masses that contain material from a small segment of 
the spiral arm (Boley et al. 2010). 

Because spiral arms are the most likely regions for fragmentation, one must also consider 
enrichment of solids at birth. Solids are expected to move relative to the gas due to the lack of 
pressure support (Weidenschilling 1977). Depending on the direction of the pressure gradient, 
this disparity between gas and particle forces can lead to either a tail or headwind for a given 
particle. The wind exchanges angular momentum with the solid and causes it to migrate. 
Haghighipour & Boss (2003) showed that this effect will cause particles to become trapped 
in local pressure maxima, and Rice et al. (2004, 2006) used SPH calculations with solid 
particles to demonstrate that spiral arms in a gravitationally unstable disk are ideal places 
for trapping solids. These results suggest that enrichment at birth could provide newly- 
formed fragments with a supernebular solids-to-gas ratio. In addition, such concentration 
of solids could provide enough material to form substantial rocky/icy cores, regardless of 
whether the fragment survives in the disk. 

In this paper, we present hydrodynamics+particle simulations to explore the heavy- 
element enrichment of fragments at birth in spatially extended, massive disks during the 
accretion phase. In Section 2, we outline the numerical methods for this study, including 
introducing our algorithm for evolving solids. We then describe the numerical experiments 
and disk initial conditions in Section 3. We present our results in Sections 4 and 5, and 
discuss the implications of some of these results in section 6. We summarize our conclusions 
with a buUeted list in Section 7. 



2. Method 

2.1. Hydrodynamics 

We use CHYMERA (Boley 2007), with the particle modifications described below, to 
run a series of hydrodynamics simulations of gravitationally unstable disks with a signif- 
icant population of rocks. CHYMERA is an Eulerian code that solves the equations of 
hydrodynamics, with self-gravity, on a regular, cylindrical grid. The star's motion is solved 
self-consistently, as described in B2009. The internal energy of the gas accounts for the 
vibrorotational states of H2 as described in Boley et al. (2007a). Because we are studying 
extended disks that have very long orbital timescales, we assume that the ortho-para ratio 
of H2 remains in equilibrium for the local temperature. For radiative cooling, CHYMERA 
can use flux-limited diffusion (Boley et al. 2006), a hybrid ray-tracing+flux-limited diffusion 
scheme (Boley et al. 2007b), or a local radiative cooling approximation (B2009). The hybrid 
scheme is the most accurate, but because it solves for radiative transport explicitly, it can 
become unstable whenever the radiative time step becomes smaller than the hydrodynamics 
step. Under most conditions, the algorithm is stabilized by limiting the amount of energy a 
cell can change for a given step. This sacrifices accuracy for speed, but appears to be reason- 
able for many situations (see Boley 2007). To avoid using limiters, a subcycling version of 
the hybrid scheme has been developed, but it is still in its testing phase and was not ready 
for use when this study began. Because we are interested in how the solids distribution 
behaves in a disk that does fragment, we choose to use the fast and stable local radiative 
cooling approximation. 

In order to suppress artificial fragmentation, we ensure that the Truelove et al. (1997) 
criterion (see also Nelson 2006) is obeyed at all times by adding a cold pressure floor to our 
equation of state, i.e., P = max ( Pineal > -Pcow)- We have chosen to set the minimum Jeans 
wavelength to be at least four radial grid cells, which sets Pcow = (4Ar)^G'p^/7r. As will be 
discussed in more detail below, the cold pressure only becomes necessary at the very center 
of clumps. 



2.2. Gas Drag Algorithm 

For modeling solids, we have augmented CHYMERA to include particles, where each 
particle represents an ensemble of solids with some radius s. The Epstein and Stokes drags 
are calculated separately, and an interpolated solution is then found. First consider the 
Epstein limit, which is appropriate whenever the Knudsen number Kn = A/2s ^ 1. Here, 
A = fimp{nR'^p)~^ is the mean-free-path of a gas particle, where p is the gas density with 



mean molecular weight /x, rrip is the proton mass, and R is the typical gas atom/molecule 
radius. The internuclear spacing of H2 is tq ~ 0.7AA, and the typical diameter for a molecule 
is ~ 3ro (Allen 1976), so we take R = 10~^ cm. This size-scale is slightly smaller than what 
has been used in other work (cf. Rice et al. 2004), but gives a cross section that is within a 
factor of a few to expected cross sections for atomic and molecular gas (~ 10^^^ cm^; Allen 
1976). 

We set the velocity difference between the gas and solids to be Sv = Vg — Vp. Following 
Paardekooper (2007), we write the change of this difference in the Epstein limit as 

at Epstein V^/ PsS 

for adiabatic sound speed Ca and internal particle density ps- The term fe is used to interpo- 
late smoothly between low and high-Mach speed flows, where the Mach speed here is defined 
relative to the gas such that M = 6v/ca and /^ = (1 + 97r5v^/{l28cl)y/^ (Kwok 1975). It 
is also convenient to define the stopping time for a particle moving through the gas as (e.g., 
Weidenschilling 1977) 

*.= ^. (2) 

pCa 

Assuming that the gas density and sound speed are constant over a time step At, equation 
(1) can be integrated analytically. Carrying out the integration, we find 



Sv 



= 26voPQ exp {-aAt) (Q^ - Sv^ exp (-2aAt)) ^ , (3) 

Epstein 



where Svq is the initial velocity difference, a = (^) ^, (3 = (128c^/(97r)) , and Q 



/3 + (/3^ + Svq) (see also Paardekooper & Mellema 2006 for an equivalent expression). 

Now consider the Stokes limit, which is appropriate for small Kn. The change in the 
velocity difference can be written as 



d6v 
~dt 



= -SaKnkdSv. (4) 

Stokes 

The factor ka depends on the Reynolds number Re = 3 (|) ^, where (e.g., Paardekooper 
2007) 

l + 0.15Re°-^®^ for Re < 500, 

kd = { 3.96 X lO-^Re^-^ for Re < 1500, (5) 

0.1 IRe otherwise. 



Each Reynolds regime can be integrated separately, yielding the following solutions: 



5v 



5t;oe-3-KnAt (e|5^p|0.687 (^^ _ g-2.061aKnAt) ^ ^^j -1-4556 ^^^ ^^ ^ ^^^^ 

6vo (1 + 7.2\6vo\^-^KnaTAtf°-^^^'^ for Re < 150( 

6vo (l + 0.99 (f ) ^^^ ^At) "^ otherwise. 



(6) 



For convenience, we have defined 6 = 0.15 ( (|) ^^ ] and T = 3.96x 10"^ f (|) 



Ca Kn 



which are constant over a time step. 



Finally, the Epstein and Stokes solutions can be coupled to find an interpolated solution 
(modified from Woitke & Helling 2003)0 



Sv 



\2 



(=*«"- .,„ 



total (3Kn)^ + 1 Epstein ' (3Kn)^ + 1 



1 



(7) 

Stokes 



Using conservation of momentum, this new velocity difference can be used to find the new 
momentum and velocity for the gas and solids. We do include the back reaction of the 
particles on the gas in this algorithm. 

We should point out that a limitation of our approach is that the analytic solutions 
for the velocity difference consider only the drag force. However, differences between forces 
on the gas compared with those on the particles due to pressure require that either these 
forces be included in the solution of the velocity difference or that the hydrodynamics time 
step remains smaller than the stopping time of the particles. If these requirements are not 
met, then the asymptotic velocity drift of the particles through the gas will be modeled 
poorly. For the particle sizes and density regimes we consider here, this should not be a 
major difficulty. We present a drift test in the Appendix that addresses this issue. 



2.3. Coupling Gas and Solids 

We use a crude particle-in-cell method for evolving the particles. At the start of the 
simulation, particles are assigned to cells randomly, weighted by the gas mass in a given 
cell. Inside the cell, each particle is given a random position, an orbital frequency that is 

^The interpolation weights used by Woitke & Helling are (3Kn+i')2 ^^^ TsKn+TF" However, they apply 
these weights directly to the force, where we solve for the velocity difference for each limit independently 
and seek a coupled solution. For this reason, we modify the weights to sum to unity. The discrepancy in the 
sum of the weights is, at worst, different by a factor of two at Kn = 1/3. Figure 21 in the Appendix shows 
that drift velocities are captured with accuracy sufficient enough for this study. 



Keplerian, and a random vertical and radial velocity as explained in Section 3.2. Whenever 
the density of solids is needed, it is calculated by summing the mass of all particles in a 
given cell and dividing by that cell's volume. As discussed below, this density is included 
when calculating the disk's self-gravity. For future studies, we would like to include density 
estimates using a smoothing kernel over nearest neighbors, but that is not implemented here. 

We integrate the particles using a leap-frog method in Cartesian coordinates. During 
a given hydrodynamic step At, the particles are integrated after the second hydrodynamics 
sourcing step, which occurs after the potential update (see Figure 2.2 in Boley 2007), using 
the following scheme: 

^^0 = ^-1/2 + -agravAt_i (8) 

^^1/2 = ^0 + -OgravAto 

^i/2 = G'drag('^^^,^^l/2; Ato) 
Xi = Xo + f i/gAto. 

Here, v and x are the velocity and position values in Cartesian coordinates. The subscript 
represents the values at the beginning of the current step. Using the previous time step with 
the current acceleration for the first update allows for better accuracy during variable time 
stepping. The velocity is first updated using the acceleration due to gravity. Next, using the 
updated velocity, a new velocity v' is found by applying the drag force using equations (3), 
(6), and (7) and using conservation of momentum (see below), all represented in equation 
(8) as Gdrag- After the particle velocities are updated, the positions are advanced using v'. 
This algorithm permits 10^ particles to be evolved with a modest performance penalty. 

For calculating Ograv, we first apply the exact force on the particle due to the star using 
Cartesian coordinates. Next, the force due to self-gravity from the disk is calculated on the 
cell faces of the cylindrical grid. These forces are then interpolated linearly, in the same 
way as the velocities (see below), to the particle's position in cyhndrical coordinates. These 
interpolated values are then transformed to Cartesian coordinates and added to the star's 
force. The drag force is calculated by transforming the ith particle's Cartesian velocities to 
cylindrical components Vr^i, fz.j, and f^^j. The initial velocity differences used in equations (3), 
(6), and (7) are then calculated by linearly interpolating the gas velocities to the position of 
the particle. These interpolations, for now, are only one dimensional. For example, regardless 
of a particle's (p or z position within a given cell, the face-centered radial velocity component 
is used for interpolating to the r position for finding the gas Vr- Likewise, regardless of a 
particle's or r position in a cell, the face-centered vertical velocity component is used for 



interpolating to the z position for finding the gas v^. When calculating the local gas orbital 
frequency Vt within the Jth radial cell, we use fi = fij + (r2j+i — r2j_i)(rj — rj)/{2^r), 
regardless of the particle's (p ot z position within the given cell, for particle radial position rj 
and cell center position rj. This interpolation scheme is coarse but sufficient for this study 
(see Appendix). Using the time step At/2, the new 6v is calculated for each component. So 
far, we only have the new velocity differences. We now use conservation of momentum to find 
the new velocity for the particle and for the gas, as affected by this particle. For example, in 
the radial direction, we have 6vrfi = i^r,g,o — i^r,p,i,o and initial total radial momentum density 
S = PgVr,gfi + Pp,iVr,p,i,o, wherc Vr,gfi IS the current gas radial velocity at the interpolated 
position, pg is the gas density for the cell (cell center value), and fr,p,i,o is the ith particle's 
radial velocity. The density for the ith particle, Pp^i, is found by dividing the particle's mass 
by the volume of the cell that contains the particle. The naughts represent values before 
the drag from the ith particle is applied. The updated 6vr is then found by using Svrfi in 
equations (3), (6), and (7). The new velocity in the radial direction for the ith particle is 
then Vr^p^i = {S — Pg6vr) / {pg + Pp^i) . The new gas momentum at the location of the particle is 
then calculated. In the radial and vertical directions, the momentum densities are updated 
on cell faces by extrapolating the change in the gas momentum density to the cell faces using 
the same interpolation weights that were used to calculate the gas velocity at the particle's 
position. This is done to maintain the staggered grid formalism used in the code. Because 
the angular momentum density is cell centered, the change in the gas's angular momentum 
is applied directly to the cell in which the particle is contained. Finally, we should make 
clear that every time the drag force is applied, the velocities for both the gas and the ith 
particle are updated. If the ith+1 particle resides within the same cell as the ith particle, 
then the drag on the ith+1 particle will be affected by the drag update for the ith particle. 
Waiting to update the gas velocities until all particles have had their drag forces taken into 
account leads to numerical instabilities when the particle density approaches the gas density. 



3. Initial Conditions 

3.1. Gas Model 

A principal objective of this study is to address whether clumps that form by disk 
instability can be enriched at birth. To proceed, we explore two disks, each of which is 
forced to become highly unstable in order to effect fragmentation in the system. The first 
initial condition (ICl) is a 0.4 Mq disk around a 1.5 Mq star. We use a softening for the 
star of 5AU as described in B2009. The ICs were created using the analytic approximations 
described in Boley & Durisen (2008). Unfortunately, this approximation uses Keplerian 



rotation, assumes a polytropic vertical temperature profile, and has sharp inner and outer 
boundaries. As a result, the analytic disk was evolved at low azimuthal resolution for about 
two orbital periods at r ~ 200 AU to allow the disk to adjust. During this evolution, the 
disk was irradiated using a temperature profile of Tjrr = To(r/AU)~^/^ + 10 K, with Tq = 600 
K. Figure 1 shows the corresponding Toomre (1964) Q and surface density S profiles. Recall 
that a disk is unstable to gravitational instabilities whenever Q = CaK/{7rGT,) < 1.7, where 
K is the epicyclic frequency and Ca is the adiabatic sound speed (e.g., Durisen et al. 2007). 
In general, instability does not mean fragmentation, as GIs are capable of reaching a state 
of self-regulation (e.g., Gammie 2001; Lodato & Rice 2004; Mejia et al. 2005), even when 
radiation transport is included (Nelson et al. 2000; Boley et al. 2006, 2007; Stamatellos 
& Whitworth 2008). However, at least for isothermal disks, fragmentation becomes very 
likely whenever Q < 1.4 (e.g. Nelson et al. 1998; Mayer et al. 2004). Figure 1 shows 
that the high irradiation temperature used for the initial model stabilizes the disk against 
fragmentation, but not against the development of gravitational instabilities (GIs). During 
the low-resolution evolution, the surface density profile has relaxed to an exponential that 
roughly follows S ~ 410exp(— r/45AU). This is not a detailed fit, as we are looking at 
a small range of disk radii, but it does give the basic sense of the surface density profile. 
Hereafter, simulations that use ICl are labeled SIMl. 

At the start of the simulation, SIMl's ICs are loaded onto a high-resolution grid, with 
r, z, (f) = 256, 64, 512 zones. The initial radial and vertical cell widths remain unchanged 
from the relaxation process, with Ar = Az = 2 AU. The disk is given a 10% random 
density perturbation, and the prefactor Tq in the irradiation temperature profile is reduced 
to 130 K. Dropping the irradiation temperature by such a large factor causes a sudden loss 
of disk stabilization and allows Q to drop below unity temporarily. The decision to change 
To dramatically is mainly to drive SIMl to instability violently; however, it does have the 
physical interpretation of sudden and prolonged shielding of the outer disk by, for example, 
the inner disk. 

The second set of ICs (IC2) uses the B2009 simulation called SIMD. We start our 
realization of the simulation about half of an orbit before fragmentation (Fig. 2). To reca- 
pitulate, the disk is 0.33 Mq and surrounds a 1 Mq star. We do not use softening for the 
interactions between the star and gas in this simulation. The grid has the same dimensions 
and physical resolution as the SIMl grid. The disk is grown by allowing mass accretion 
from an envisaged envelope at M ~ IO^^Mq yr^^ until 0.33 Mq is reached, after which the 
accretion is abruptly halted. This mass accretion rate is high, and is chosen mainly as a 
computational time consideration (see B2009 for more details). The radial density profile for 
the infalling mass is set to p ~ j--o.5 (^^ote that r~^'^ profiles were also explored in B2009). 
This shallow density profile causes mass to build up at large radii, which creates a surface 
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density enhancement in the disk's profile (Fig. 3) and causes Q to plummet. The surface 
density profile reflects that, for this particular case, disk formation occurs faster than Gls 
can activate and smooth out the density structure. During the formation of the disk and its 
subsequent evolution, an irradiation temperature of 30 K is incident on the disk everywhere. 
We include this irradiation for the simulation presented here, which we call SIM2. 



3.2. Solids and Opacity Assumptions 

We assume that half of the solid mass is contained in large particles that do not con- 
tribute to the opacity and that half goes into a perfectly gas-entrained dust distribution that 
is solely responsible for the opacity. We use the D'Alessio et al. (2001) opacity tables, and 
only consider the Rosseland mean for simplicity. The grain size distribution is assumed to 
be a power law, with dn/ds oc s~^'^ and a minimum grain size of 0.005/xm. The maximum 
grain size in the distribution can be varied, which can have a strong effect on the cooling 
(e.g., Cai et al. 2006). SIMl is run with a maximum grain size of 1 mm and with 1 /im, 
while SIM2 is only run with 1 /im. Hereafter, when we use the term "1 mm opacities," we 
are referring to a dust grain distribution with a maximum s = 1 mm. The total solids- 
to-gas mass ratio is assumed to be 0.02 (including ices), with 0.01 in small grains, which 
sets the opacity, and 0.01 in rock-size solids composed of silicates and ices. We refer to the 
solids in this manuscript simply as "rocks" because most simulations are run with lOcm-size 
particles, but rocks can also refer to the solids in the Ikm-size run. We choose this simple 
division between large solids and grains for multiple reasons, which are detailed at the end 
of this subsection. Note that SIMD in B2009 was assumed to have solar metallicity, so our 
assumption here reduces the original dust opacity for SIM2 by a factor of two. 

Particles are distributed according to a density-weighted Monte Carlo sampling of the 
disk. The initial particle surface densities are plotted in Figures 1 and 3 for SlMl and S1M2, 
respectively, and show good agreement with the gas surface density. We assume that each 
particle in the simulation represents an ensemble of solids. The particles in any given run are 
assumed to have a single size population of either 10 cm or 1 km in radius and an internal 
density of 3 g cm^'^. For the initial velocity distribution, the particles are placed on Keplerian 
orbits for the star. The disk potential in the initial orbital frequency was not included because 
we intended to focus on just 10-cm size particles, which respond strongly to gas drag. The 1 
km-size run was originally intended to be presented in a complementary follow-up study, but 
for providing a comparison with the low-drag limit, it was later decided to run SIM2 with 
km-size particles. To make sure that any difference we find can be attributed to the change 
in solid size only, we also put the 1 km-size particles on Keplerian orbits. In SIMl, the radial 
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and vertical velocities were set to zero and given a random perturbation with a magnitude 
no greater than 0.01 of the Keplerian speed. The perturbation was included in an attempt to 
avoid all particles at a given disk radius from passing through the midplane simultaneously. 
In SIMl, we found the 0.01 perturbation to be too weak, so the perturbation was increased 
to 0.1 of the Keplerian speed for SIM2. The larger perturbation kept strong caustics from 
forming during particle midplane passages. However, the fast midplane settling times (~ 1 
orbit) marginalize this detail for the 10 cm-size solids. 

For this study, we do not change the fraction of solids that go into dust grains. Although 
our 50-50 choice is arbitrary, the effects of decreasing or increasing the amount of dust in the 
system on GI amplitudes have already been addressed in previous research (see Cai et al. 
2006, 2008; Boley & Durisen 2008). As more solids remain in a dusty phase, the opacity of 
the disk will increase, which will tend to stabilize the disk against fragmentation. Likewise, 
as the opacity is lowered due to more dust going into larger solids, the amplitude of GIs 
can increase due to shorter cooling times, as long as the emissivity does not approach zero. 
Another consideration is whether there is a qualitative change in the disk's behavior for 
different large solids-to-gas ratios, which could occur for different disk metallicities or for 
different grains/large solids fractions. The 50-50 choice seems to be a reasonable starting 
point, as it is not at either extreme, and, as we will show, does demonstrate enrichment and 
differentiation at birth. To begin to understand how the large solids affect the large-scale 
structure of the gaseous disk, we run one simulation without large solids, while keeping the 
opacity the same. 

We explore two particle sizes, but not a distribution of particle sizes. There are several 
reasons for this choice. First, we seek to run as few supercomputer simulations as necessary 
to demonstrate the enrichment of fragments at birth. We feel this can be done well by 
picking a solid size that is likely to maximize enrichment and differentiation at birth, which 
is represented by 10 cm rocks (see the appendix). Particles that are much smaller will 
remain entrained with the gas, while particles much larger will be unaffected by gas drag. 
The former does not require additional simulations, as this has been the assumption in 
most GI hydrodynamics simulations to date. The latter is demonstrated using the 1 km 
solid size simulation, but it is not expected to lead to significant enrichment (e.g.. Helled 
& Bodenheimer 2010). Instead of exploring a single population, one might argue that a 
distribution of particle sizes would be more realistic. While this is correct qualitatively, any 
distribution that we could have chosen would be as equally arbitrary as a single solid size 
inasmuch as the actual distribution of solid sizes is unknown, particularly during the early 
stages of disk evolution explored here. Even during later stages, planetesimal growth is a 
matter of ongoing research (e.g., Johansen et al. 2007; Morbidelli et al. 2009). Consider, for 
example, the results of the Brauer et al. (2008) study. They find that the particle size of 
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solids can become strongly peaked between 1 and 10 cm at 100 AU for the conditions they 
explore, similar to what is explored in this manuscript, but only if they ignore collisional 
destruction as well as inward drift due to gas drag. If they include destruction and drift, 
they find very little dust growth, which is inconsistent with observations. For example. 
Greaves et al. (2008) and Ricci et al. (2010) both find evidence for particles with sizes of 
a few cm or larger in the extended regions of disks. Even if we use these observations to 
infer a power law distribution for particle sizes up to 10 cm (e.g., Draine 2006), we are still 
are left without knowing the distribution of solids larger than 10 cm, due to observational 
limits, or whether the distribution appropriate for disks with ages of a few Myr is also 
appropriate for very young disks. Finally, there is a danger in assuming some distribution 
of solids without having a model for grain growth and destruction that can be evolved self- 
consistently. As solid sizes differentiate due to disparate drag timescales, the local population 
of solids will also evolve. The link between the spatial and size evolution of the solids will 
alter the differentiation process itself, as well as local opacities. This problem is by no means 
removed by considering a single particle size, as a local concentration of rocks will likely 
change the local dust opacity, but a single size does allow as clean of a numerical experiment 
as possible. Complexity can be added in later studies after the simple case is understood, 
which we present in detail here. 



3.3. Simulation Nomenclature 

Overall, four variations of the simulations are run. As is described in the next sec- 
tion, the SIMl disk is evolved with 1mm opacities, with lyum opacities, and the same 1/im 
opacities, but without rocks. We refer to these simulations as SIMlmm, SIMlmu, and 
SIM ImuNo Rocks, respectively. SIMlmu and SIMlmuNoRocks are both branched from the 
SIMlmm disk after the first 620 yr, as described below. When referring generally to char- 
acteristics of the disk, we only write SIM. Two realizations of SIM2 are presented in this 
study, one using 10 cm-size particles (simply SIM2) and 1 km-size particles (SIM2km), both 
of which are run with l/im opacities. 



4. Results 

4.1. SIMl: Gas Evolution 

Due to the short cooling times, the sudden change from Tq = 600 K to 130 K causes the 
disk to contract quickly, and after about 620 yr of evolution, the SIMl disk has contracted 
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enough that it becomes possible to evolve the same disk on a grid with half the initial physical 
size in r and z. Using direct injection, i.e., no interpolation, the disk is loaded onto a higher 
resolution grid with Ar = I\z = 1 AU, while keeping the number of radial, azimuthal, and 
vertical grid cells the same. Figure 4 shows the Q and surface density profiles shortly after 
being loaded onto the high physical resolution grid. The surface density profile has been 
altered significantly after reducing the irradiation temperature. The choppy Q profile is due 
to direct injection, which adds noise to the epicyclic frequency. According to Figure 4, the 
disk is strongly unstable. 

SIMlmm is only run for 1260 yr. The surface density evolution is shown in Figure 3. 
Although very flocculent structure develops, including knot formation, clumps do not form. 
If, however, the maximum grain size is reduced to 1 /xm at 620 yr, which marks the start 
of SIMlmu, the disk fragments into many clumps. At the time of the opacity switch, dense 
spiral arms have not yet fully developed. Figure 6 shows the evolution of SIMlmu from 620 
yr to 1870 yr. The effect of opacity on the spiral structure is illustrated by the p-T phase 
diagram in Figure 7. For the same temperature, SIMlmu reaches much higher densities than 
SIMlmm. Figure 8 shows the D'Alessio Rosseland mean opacity for a maximum grain size of 
1 /xm and 1 mm. For T < 100 K, the 1 mm opacity is considerably more opaque, which does 
not allow the spiral arms in SIMlmm to cool fast enough for fragmentation to proceed. One 
might worry, instead, that the cold pressure (see section 2.1) employed in these simulations 
could be affecting fragmentation. The thick line in Figure 7 shows where the cold pressure 
becomes important, assuming an equivalent thermodynamic temperature for the pressure. 
This threshold is only met in clump cores, so this cold pressure cannot be responsible for 
the lack of fragmentation in the 1 mm case. It is possible that, after significantly longer 
evolution, SIMlmm would fragment, too. We chose not to investigate this possibility for 
two reasons: (1) We wanted to focus our computing resources on disks that we know will 
fragment; and (2) even this short simulation demonstrates that the grain size in the outer 
disk can strongly affect clump formation. 

Finally, before focusing on the evolution of the rocks in these simulations, we compare 
disk structures between SIMlmu and SIMlmuNoRocks (Fig. 9). Although the disks are 
qualitatively similar, differences have developed. This will be discussed in more detail in the 
next section. 



4.2. SIMl: Rock Evolution 

By the time spiral arms in the gas form, most rocks have settled into a very thin subdisk, 
contained within one vertical grid cell. This is consistent with expectations: As discussed in 
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the Appendix, the stopping time for 10 cm rocks is ~ 10-100 yr for the gas densities in these 
disks. Most particles will also make their first midplane passage after 1/4 of an orbit due to 
the ICs, i.e., about 200 yr at 100 AU. As seen in Figure 6, the spiral arms remain weak until 
about 780 yr, giving enough time for particles to settle before spiral waves fully develop. 

Figure 10 shows the face-on evolution of the gas and particles in SIMlmu. As rocks 
encounter gas over-densities in spiral arms, they become trapped and concentrated, which 
is consistent with the results reported by Rice et al. (2004). By the end of the simulation, 
most rocks are associated with spiral arms or clumps. To determine whether there is an 
enhancement of solids in these structures, we plot in Figure 11 the ratio of rock surface 
density to gas surface density for four snapshots. The data are shown using the cylindrical 
coordinates of the grid, where the abscissa is the r coordinate and the ordinate is the 
coordinate, to marginalize the chance that a high peak in the ratio is due to interpolation 
error. In about one dynamical time, the solids become filtered and caught in the spiral arms, 
with some regions exceeding a tenth of the gas surface density. The right panel in Figure 11 
shows the rock-to-gas ratio for gas and particles in only the midplane cells, which accounts 
for most of the rock mass. Thin, wispy rock arms dominate the midplane density in some 
regions. 

In addition to concentrating the rocks into arms, the gas-rock coupling in this highly 
perturbed disk stops rapid inward drift. Rocks still move inward, but they are locked with 
the progression of spiral arms. Figure 12 (left panel) shows the azimuthally averaged surface 
density evolution for the rocks. For the first ~ 600 yr of evolution, particles outside r ~ 70 
AU move inward relative to the gas, while particles inside 70 AU move outward and collect 
in an annular surface density maximum. After ~ 600 yr, spiral arms develop in the 70 
AU surface density enhancement, and undermine simple inward drift due to nontrivial gas 
dynamics and the evolving pressure maxima. Note that rocks move inward and outward 
over tens of AU. The spikes in the rock surface density correspond to gaseous spiral arms 
and/or clumps. As shown in Figure 12 (right panel), the surface density enhancements in 
rocks coincide with maxima in the azimuthally averaged gas surface density. 

Occasionally, features that do not have an obvious gas surface density counterpart form. 
The most obvious type of feature is a rock bridge between gaseous spiral arms. These 
structures are created when a gas spiral arm dissipates, but the thin, dense rock arm does 
not. As the spiral arm dissipates, it does so by diffusing away from the pressure maximum. 
This would allow for particles that are highly concentrated and not perfectly entrained with 
the gas to remain concentrated. Bridges mostly have irregular shapes, but rare, straight 
bridges do form after interacting arms create a temporary leading rock arm that is then 
sheared into a trailing arm. Knotty structure in the rocks can also form due to buckling of 
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the rock arm after incomplete dissolution of the gaseous arm (see x,y ^ 50,50 AU in the 
1410 yr panel of Fig. 10). 

Finally, the strong clumping and density enhancements in the rock distribution seem to 
be causing noticeable changes in the evolution of the gas. In particular, as seen in Figure 
9, SIMlmu shows additional clumping and fragmentation compared with SIMlmuNoRocks. 
Unfortunately, it is difficult to distinguish whether the differences between the simulations 
are brought about by the chaotic nature of GIs in the nonlinear regime, or whether they are 
directly related to the solids. This is particularly difficult because variations of a percent 
could be amplified over a dynamical time for either reason. For example, by removing the 
solids from the simulation, we instantaneously decrease the mass of the disk by one percent. 
On the other hand, significant concentration in spiral arms could create small, but important, 
gravitational potential enhancements, aiding fragmentation. Additional effects, such as the 
momentum exchange between the gas and the rocks, may have some effect as well (e.g., 
Johansen et al. 2007). Moreover, a greater fraction of large solids will likely exacerbate the 
differences, if they can indeed be attributed to solid-gas interactions. This result is a topic 
by itself, and needs to be addressed in a separate study. 



4.3. SIM2: Gas and Rock Evolution 

SIM2 allows us to explore fragmentation and concentration of solids in a disk that is 
driven toward instability by mass loading. It is not as violently unstable as SIMl, and has 
much more ordered spiral arms. In addition, visible, weak spiral structure is already present 
at the start of the evolution presented here. Because the particles are distributed according 
to gas mass, the rocks have not fully settled to the midplane by the time dense spiral arms 
develop, unlike SIMl. This allows for the possibility that the spiral waves might significantly 
hinder midplane rock settling. 

Within about 1/2 orbit at r ~ 100 AU, an arm becomes very dense and fragments into 
a clump. The simulation is evolved for an additional 1/2 orbit after fragmentation. At the 
time the simulation is stopped, a second arm is showing signs of clump formation, but we do 
not investigate it further here. The surface density evolution of the disk is given in Figure 
13. As done for SIMl, the same surface density snapshots are shown in Figure 14, but with 
the particles superimposed. The rocks respond quickly to the density perturbations, and the 
spiral arms become enhanced despite the short duration of the simulation, as highlighted in 
Figure 15. Such rapid concentration indicates that enhancement of solids will occur as soon 
as spiral arms form. Even with large-scale spiral arms, significant settling to the midplane 
takes place. Within one orbit, the midplane has a very high solids concentration, with 
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> 95% of rocks mass contained in the midplane cells. In some spiral arms, the rock density 
dominates over the gas density. 

In Figure 16, we show the azimuthally averaged rock surface density profile for several 
snapshots. Particles move inward, but are thrown back out to larger radii as the spiral arms 
increase in amplitude. As seen in SIMl, rocks are concentrated in azimuthally averaged 
density maxima. In contrast to the 10 cm-size particles, the solids in SIM2km do not collect 
in the spiral gaseous spiral arms (Fig. 17). Spiral arms do form in the solids, but because the 
1 km-size solids respond to gas drag slowly, the particles do not collect in the high-density 
gas regions. This, in turn, allows for super nebular solids-to-gas ratios outside of the gaseous 
spiral arms (Fig. 18). The enhancement is partly due to the much broader solid spiral arms, 
causing enhancement in the low-density regions immediately before and after a gaseous arm, 
and partly due to misalignments between the gaseous and solid arms. 



5. Clumps 

One of the principal science questions that we want to address in this study is whether 
clumps that form by disk instability can be enriched at birth. As shown in the previous 
sections, spiral arms sweep up rocks rapidly, and whenever fragmentation does occur, clumps 
are born in these rock-rich environments. In this section, we address the magnitude of this 
enrichment and the consequences it has for gas giant composition and core formation. 

Consider first the clump that forms in SIM2 (SIM2:C1) because it has a relatively clean 
formation environment, i.e., it forms from one spiral arm and does not suffer mergers or 
collisions with other arms or clumps. Figure 19 shows the rock-to-gas density ratio /ntoG for 
the full surface density (left) and for just the midplane cells (right). The center of the clump 
{(f) ~ 5.3 rad) shows a large rock enhancement, as do the spiral arms/wakes that continue 
to supply solids and gas to the clump. The dearth of rocks at larger azimuth (prograde 
direction; see Fig. 14) is not completely understood, but seems to be related to formation of 
the nearby structure. 

In order to quantify the rock enrichment, we need an estimate for the clump mass. 
Figure 20 plots cell gas volume density and the gas temperature as a function of radius from 
the densest cell in the clump. The minor peaks in the profiles owe to the highly fiattened 
clump structure, yielding poor vertical resolution after collapse. Based on these profiles, we 
consider a given cell's mass to belong to the clump if its gas density p ^ 9 x 10~^^ g cm~^"^ and 
its temperature T > 34 K. Although subjective, these cutoffs provide a reasonable definition 
for the clump, and are chosen to be low, but still able to distinguish the clump from spiral 
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arm material. The clump masses vary by a few percent when varying the cutoff temperature 
by a few Kelvin or the density by ~ 10%. Figure 21 shows a surface density+rocks snapshot 
of the clump at 1300 yr (right) and just before fragmentation at 960 yr (left). If a particle is 
located in a cell that belongs to a clump, as defined above, we consider the particle to also 
belong to that clump. This criterion will not work in general, but for the 10 cm particles 
followed here, the criterion is reasonable. Particles that belong to the clump are plotted 
using blue crosses, and seem to delineate the clump's extent well. These particles are also 
shown on the 960 yr panel, which demonstrates that the clump forms from a 20-30 AU 
section of the arm, consistent with Durisen et al. (2008) and Boley et al. (2010), with some 
gas and rock mass growth after its initial formation. The clump at this stage of its formation 
is deformed, and although the solids are not concentrated near the geometric center of the 
clump, they are concentrated at the clump's density maximum. 

Table 1 gives the gas and rock mass of SIM2:C1 at 1300, 1450, and 1620 yr. We also 
show mass estimates for gas with p > pmax/2 and for p > 0.9pmax- These additional cutoffs 
allow us to explore the rock distribution in the clump. Over 320 yr, the clump grows from 7.1 
to 8.5 Mj. The rock mass also increases from 35 to 42 M^, keeping f-RxoG ~ 0.015 throughout 
the clump's evolution. When taking into account the entrained dust distribution, the overall 
enrichment of the clump /enr = {fmoG + 0.01)/0.02 ~ 1.3. Despite forming from the spiral 
arm, the enrichment of the entire clump is fairly modest, about 30% over the nebula's average 
value. We can understand this result as follows: Although the center of the spiral arms are 
highly rock-enriched, the arm material immediately surrounding the midplane of the spiral 
arm is rock depleted. Collapse mixes this material back together, resulting in a smaller 
enrichment than what one might expect from the densest regions of the spiral arm alone. 
Mass growth does not appear to change this ratio significantly. This is likely due to growth 
being mediated by filaments (see Fig. 19), which are rock-enriched. 

Even though the clump as a whole has modest enrichment, the rocks are highly con- 
centrated near the center. At half the peak density, the clump is enriched by a factor of a 
few in the rocks and about two in total solids. Near the center of the clump, the rock-to-gas 
ratio is ten times the average nebula value, with 30 to 38 M^ of rock concentrated at the 
center. Even if the entire clump is eventually disrupted, extreme rock concentration could 
lead to the formation of early cores and planetoids in the outer disk. 

In contrast to SIM2:C1, SIM2km:Cl is depleted in rocks, as km-size solids are not as 
easily captured as the 10 cm-size solids. Table 1 shows that there is no additional concen- 
tration of solids toward the center, so unlike the 10 cm case, no differentiation has taken 
place. Overall, the clump is subnebular in solid composition, is less massive than its 10 cm 
counterpart, and has a lower central density. This provides additional evidence that the 



evolution of the solids affect the evolution of the gaseous disk, even for features as large as 
fragments. The change is not as drastic as that seen between SIMlmu and SIMlmuNoRocks, 
but SIM2 is not as flocculent as SIMl, either. 

The low level of enrichment cannot be attributed to lack of interactions between solids 
and the gaseous fragment. Figure 22 traces the solids associated with SIM2km:Cl at about 
1720 yr, determined the same way as for SIM2:C1, and traces the particles backwards and 
forwards in the clump's orbit. In the earliest snapshot shown (1290 yr), roughly 1 Mq of 
solids have yet to enter the planet's Hill sphere, while for the oldest snapshot (2190 yr), 
about 0.7 M® of the solids escape have passed out of the new, larger Hill sphere. This does 
show that some planetesimals are being captured into the orbit of the protoplanet, but it 
does not lead to enrichment. In the latter snapshot, the clump has grown to roughly 12 Mj, 
but only has 11 Mq of large solids, giving a rock-to-gas ratio of 0.0029, which makes the 
protoplanet more anemic than seen for the snapshot in Table 1. The capturing of solids for 
this clump appears to be due to the capture of only the lowest relative velocity solids, aided 
by a gradual deepening of the protoplanet's potential well due to gas accretion. 

Unlike the SIM2 simulations, SIMlmu shows prodigious clump formation with very 
flocculent spiral structure. Clumps merge, scatter, and interact with material arms. These 
interactions sometimes lead to mass growth, while at other times, the clumps become sheared 
away. For this discussion, we consider three clumps that have different but illustrative 
histories. We denote them as SIMlmu:Cl, which becomes tidally disrupted, SIMlmu:C2, 
which grows to a brown dwarf mass scale, and SIMlmu:C3, which stays a high-mass gas 
giant (Tables 2 and 3). The same density cutoff as used for SIM2:C1 provides a reasonable 
differentiation between clump and arm material, but a slightly higher temperature threshold 
(T > 40 K) is chosen. 

SIMlmu:Cl forms between 940 and 1100 y, but only survives about 1/2 of an orbit before 
it migrates inward and becomes tidally disrupted. Its density and temperature profiles are 
not shown, but its mass is given in Table 2. Even though the clump becomes disrupted, 
it is worth noting that the central concentration in SIMlmu:Cl is 27 times the nebula's 
average rock-to-gas ratio, placing 38 M^ of rock at the fragment's center. Unfortunately, 
we do not have the resolution to follow the formation of these rocks into a self-bound core, 
but the results of Helled & Schubert (2008) do suggest that formation of such a core is quite 
possible (see also Boss 1998). In particular, our clumps start differentiated, whereas Helled 
& Schubert assume homogenous composition, so the actual timescale for core formation may 
be shorter than found in the Helled & Schubert simulations. 

The second clump to form in SIMlmu, SIMlmu:C2, survives until the end of the simu- 
lation. Its temperature and density profiles for several snapshots are given in Figure 23. The 
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clump shows prodigious mass growth from 7 to 32 Mj, mainly in large bursts due to mergers 
with other clumps and from collisions with material arms. Even with its violent history and 
rapid mass growth, its total enrichment factor remains /cnr ~ 1.4 times the nebula's average 
value. This is likely due to its amalgamation with other rock-enriched structures, keeping 
/moG roughly constant. As with SIMlmu:Cl, the rocks are highly concentrated and could 
form a large core. How such initial differentiation would affect early brown dwarf evolution 
observationally is unclear. 

Finally, SIMlmu:C3 also survives until the end of the simulation, with a mass that 
hovers around 10 Mj. Unlike SIMlmu:C3, it does not merge with another clump or pass 
through a material arm. Its /enr ~ 1-4, and it has a highly concentrated rock core. 



6. Discussion 

Fragmentation at large disk radii should occur for some systems during the mass loading 
phase of disk evolution, i.e., while there is still significant infall onto the disk. Formation 
by disk instability need not exclude formation by core accretion, and both formation chan- 
nels could operate in the same system. Although the two main formation mechanisms may 
be discernible through statistical trends, as discussed in the introduction, orbital migration 
and/or scattering will make it difficult to attribute a formation mechanism to an individ- 
ual planet. Additional observational signatures are therefore desirable, and chemical and 
structural (cores) variations may allow for such differentiation. Building on the work of, 
e.g., Haghighipour & Boss (2003), Rice et al. (2004, 2006), Johansen et al. (2007), and 
Paardekooper (2007), we have included solids in our three-dimensional hydrodynamics code 
to investigate how the distribution of solids, particularly within clumps, is affected by grav- 
itational instabilities. 



6.1. Spiral Arm-Solids Coupling 

GIs form dense spiral structures that are very efficient at concentrating solids with sizes 
~ 10 cm, with surface density enhancements of solids 30 times larger than the average nebula 
value, in agreement with the Rice et al. (2004) study. Near the midplane, the mass density 
can become dominated by 10 cm rocks. Super rock-to-gas ratios even occur in SIM2km due 
to the spiral arms of solids being offset and broader than the gaseous arms. We also find that 
solids do not simply drift inward; instead, they follow the evolution of the gas spiral structure, 
including outward migration. Simulations should be extended over much longer periods than 
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studied here so that the asymptotic behavior of GIs and sohds can be addressed, but these 
results do show that sohds need not be lost by drag and that large particles could be kept 
at large radii while GIs operate. In addition, the results of SIMlmuNoRocks suggest that 
solids can affect the gas evolution. As discussed earlier, this result requires further study, 
as numerical chaos could be contributing to some of the differences that we see between 
SIMlmuNoRocks and SIMlmu. However, SIM2km also showed differences in the evolutions 
of clumps SIM2:C1 and SIM2km:Cl, indicating that more than just numerical chaos is at 
work. 



6.2. Clump Enrichment 

Clump formation takes place exactly where rock-size solids are concentrated, i.e., the 
spiral arms. Although the midplane regions of these arms can be highly enriched, disk 
instability only leads to modest enrichment within the fragments, with a total enrichment 
of tens of percent over the nebula's average value. This value may change for different 
rock distributions, but we argue that these results provide an upper limit on the amount of 
enrichment in fragments, as spiral arms are very efficient at sweeping up 10 cm rocks for the 
conditions explored here. This is emphasized by SIM2km, which actually shows a depletion 
in solids relative to the gas mass. However, Figure 18 also shows that there are supernebular 
solids near the clump for the 1 km conditions. If solid evolution is considered, then some of 
the depletion seen in SIM2km may be marginalized. Using the limits that the 10 cm and 1 
km simulations provide, we can provide an upper limit to the amount of enrichment at birth, 
assuming rock concentration does not become even more efficient for larger metallicities or 
for different opacity/large solids fractions. If we assume that most of the solids grow to 10 
cm sizes and that the rock enrichment seen in these simulations is doubled as a consequence, 
we would expect to have a typical fntoG ~ 0.03, giving a total enrichment fcm ~ 2. This 
suggests that a disk instability planet, although enriched, will have at birth no more than 
twice the amount of metals compared with its host star, taking our results at face value. If 
a substantial amount of the protoplanet's envelope is tidally stripped, the final planet could 
become more enriched than estimated from the initial clump enrichment (see section 6.4). 

Enrichment at birth provides a base level for other enrichment mechanisms, such as 
planetesimal capture. However, Helled & Bodenheimer (2010) found that enrichment by 
planetesimals with sizes > 1 km is negligible at radii of r ~ 70 AU for the conditions 
they studied. Because fragmentation is expected to be much more likely at large radii 
(e.g., Stamatellos et al. 2007; Stamatellos & Whitworth 2008; Boley 2009; Rafikov 2009; 
Clarke 2009), we expect enrichment at birth to be the dominant enrichment mechanism. 
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Nevertheless, if rapid migration occurs before dissociative collapse or if fragmentation occurs, 
however rarely, at 10s of AU, then planetesimal capture could also enhance disk instability 
objects. 



6.3. Opacity and Grain Growth 

Grain growth alters the opacity of the disk (D'Alessio et al. 2001) and changes the 
cooling rates, which in turn affects the behavior of Gls (Cai et al. 2006, 2008). Our results 
are consistent with this picture, where fragmentation is suppressed in SIMlmm due to the 
higher opacities. The problem is, unfortunately, complicated by the competing effects of grain 
growth, which suppresses fragmentation at large radii, and settling, which would reduce the 
opacity for a large volume of gas. This emphasizes that, to address this problem properly, 
grain growth models need to be included with the dynamical evolution of solids in order to 
model fragmentation of gaseous disks correctly. Concerning grain growth, we have assumed 
for this study that 10 cm or 1 km particles are present by the time the disk fragments. 
Whether this can be achieved is unclear, but the observations of Greaves et al. (2008) do 
suggest that 10 cm particles are present in massive, extended disks. 



6.4. Clump Cores and Fates 

We have assumed that the clump masses calculated here are comparable to the final gas 
giant /brown dwarf masses. Based on the arguments presented in Boley et al. (2010), we do 
not expect rapid mass growth after a clump undergoes its second collapse due to dissociation 
of H2 (the fragmentation process being the first collapse), which should occur after a few 
thousand years for these clump masses (see also Helled & Schubert 2008). Because the 
centrifugal radii of gas fragments are expected to be a few tenths to an AU at disk radii ~ 100 
AU, a substantial disk will likely form as the fragment collapses. Subsequent mass growth 
will be mediated by subdisk evolution, which has recently been demonstrated in simulations 
by Vorobyov & Basu (2010). Moreover, the clump masses calculated here represent objects 
that will become either gas giants or brown dwarfs with surrounding disks, where the disks 
themselves may be massive. Global instabilities in the clumps, e.g., the bar instability 
(Durisen et al. 1986), could also redistribute angular momentum before dissociative collapse, 
ejecting tens of percent of the clump's mass. If convection is unable to redistribute the 
centrally concentrated rocks before collapse, then such instabilities would preferentially eject 
the solid-poor material. However, even this scenario would only increase the enrichment of 
the clump by tens of percent. Only if most of the clump's initial mass could be stripped. 
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e.g., due to tides, would the enrichment of the clump match that of Jupiter's {fcnr up to 6; 
see, e.g., Baraffe et al. 2010). 

Tidal disruption of clumps with dense rock concentrations is seen in these simulations. 
For example, SIMlmu:Cl is an 8 Mj clump with 38 M^ of solids at its center, and becomes 
sheared away within half an orbit. Unfortunately, we do not have the resolution to follow 
core formation, so investigating whether these solid concentrations can separate from the 
gas (see, e.g.. Rice et al. 2006) and form large cores will need to be addressed in a future 
study. However, it is tempting to speculate whether a third formation channel for gas giants, 
which we call "core assist plus gas capture," becomes available if cores do form in clumps 
that eventually become disrupted. In this mechanism, gas fragments behave as factories 
for the formation of rocky/icy cores. After fragments become tidally disrupted, their cores 
could be liberated into the disk and begin to capture gas slowly, possibly even in the later 
and less active phase of disk evolution. This is similar to other hybrid gas giant formation 
mechanisms where hydrodynamical instabilities accelerate core formation by concentrating 
solids (e.g., Durisen et al. 2005, Klahr and Bodenheimer 2006). We also acknowledge that 
core formation through differentiation at birth may be subject to complicated equation of 
state effects as the pressure and density increase as the clump contracts. 

As we have shown here with the 1 mm opacity simulation, fragmentation is very sensitive 
to opacity, and dense knots or short-lived clumps could be more common than long-lived 
fragments (see, e.g., Vorobyov & Basu 2010). In this case, core assist plus gas capture may 
very well represent another channel for in situ formation of gas giants on wide orbits, if 
the gravitational torques and gas interactions keep the cores at large radii. If, instead, the 
cores migrate into the inner nebula as they capture gas, they will produce inner gas giants. 
Interestingly, we note that HR149026b has presented a puzzle for planet formation due to 
its extremely large core estimate (~ TOMq; Sato et al. 2005). If SIMlmu:C3 were to be 
sheared away, its concentration of TOMq of solids in its central regions could give rise to the 
formation of such a planet by core assist plus gas capture. 



6.5. Formation Channels For Substellar Companions 

Our results suggest that there may be three formation channels for gas giants planets 
(core accretion plus gas capture, direct formation by disk instability, and core assist plus gas 
capture), which could give rise to three populations of gas giants. If this is indeed the case, it 
is crucial to understand how the formation channels could be distinguished observationally. 
For this discussion, let us refer to a population I gas giant (Pop I GG) as one that forms via 
traditional core accretion plus gas capture, population II GGs as disk instability gas giants. 



23 



and pop III GGs as core assist plus gas capture. 

If pop I GGs are typically as enriched as Jupiter and Saturn, then a metallicity difference 
between pop II and pop I GGs may indeed be observable. Unfortunately, the typical outcome 
of core accretion has not yet been established, as core accretion simulations follow planet 
formation up until the runaway gas accretion phase is established, but not the end of accretion 
and subsequent evolution (e.g., see recent work by Movshovitz et al. 2010). As a result, the 
final metallicity of pop I GGs is unclear. To illustrate this point, take the enrichment of 
heavy elements in Jupiter to be about five times the solar value. If Jupiter had continued 
to accrete gas with nearly solar composition, it would need to grow to about 4 Mj before 
its enrichment would drop to twice solar and becoming indistinguishable, using metallicity, 
from pop II GGs. Looking for trends in the frequency of GGs as functions of semi-major 
axis, host star metallicity, and time, as discussed in the introduction, may still be the best 
way to distinguish pop I and II GGs. 

Pop III GGs suffer the same uncertainty in final metallicity as do pop I GGs. In the 
simulations presented here, we unfortunately cannot compute the formation of a solid core, 
let alone the subsequent runaway gas capture onto that core, and cannot state what the final 
enrichment of such an object will be. We do note however that, with the large concentration 
of solids seen in gas clumps. Pop III GGs may be very similar to Pop I GGs, but have very 
large cores for their mass (e.g., HR149026b). Whether Pop II or III objects can form in a 
given disk will likely depend on the protostellar core's angular momentum. The fraction of 
protostars for which Pop II and III GG formation will be important is uncertain, but we 
remind the reader, as also noted by Stamatellos & Whitworth (2007), that the distribution 
of angular momentum within dark cores probably does permit the formation of extended 
disks (Goodman et al. 1993). 



7. Conclusions 

We have investigated the coupling between rocks and gas in fragmenting, gravitationally 
unstable disks. We have explored the behavior of 10 cm-size and 1 km-size particles in this 
study. The former corresponds to sizes that are expected to be swept up efficiently by spiral 
arms, while the latter is largely decoupled from gas drag. Because the small-grain limit is 
simply entrained with the gas, these simulations explore the two nontrivial regimes for the 
behavior of solids in gravitationally unstable disks. Our simulations have many implications 
for disk fragmentation and planet formation, which we summarize in bullet form. 

• Spiral arms concentrate rock-size material into the regions that are most likely to 
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fragment. However, this aerodynamic capturing of solids leads to, at best, a total 
enrichment of a clump by a factor of two. A distribution of solids that does not contain 
a large amount of rock-size particles will reduce this effect, without considering solid 
evolution, as demonstrated by the 1 km-size run. 

Subsequent mass growth does not dilute the birth enrichment in these simulations 
for the 10 cm-size solids because clumps accrete gas from filaments, which are rock 
enriched, or grow by collisions with spiral arms and/or other clumps, which are also 
rock enriched. 

If most of the solid mass is in 1 km-size solids or in sizes that are equally decoupled 
from the gas, clumps can become anemic in heavy elements. 

Rock densities can become comparable to or exceed local gas densities in some struc- 
tures. Even for the 1 km-size case, regions of supernebular solids-to-gas occur due to 
differences between the spiral structures in the solids and the gas. 

The rock distribution appears to affect the behavior of gravitational instabilities and 
fragmentation. The influence of the solids on the gaseous disk, particularly fragmen- 
tation, may mainly be due to gravitational potential changes. Additional simulations 
with distributions of solids need to be explored, but must be coupled with a solid size 
evolution model. 

Gravitational instabilities interfere with the simple inward drift of solids. Instead, rocks 
that would normally have a very high inward velocity are entrained with the evolution 
of the spiral arms, and are even transported out to large radii. 

The fragmentation process is affected by the grain size distribution, i.e., the opacity. 
For the temperatures in the outer disk, larger grain sizes can stop dense condensations 
from forming long-lived clumps. 

Although the total enrichment is modest, solids are highly concentrated near the center 
of the clump as part of the fragmentation process. The solids-to-gas ratio can be as 
high as a few tenths in these regions, with tens of Earth masses of rock available for 
immediate core formation, i.e., the material does not need to settle from the entire 
clump. This could significantly decrease the core formation timescales calculated by 
Helled & Schubert. (2008). Overall, the results suggest that disk instability planets 
can have massive cores. 

Fragments can become sheared away by tides, but because rocks are concentrated near 
the center of the clump, these fragments may still form rocky cores before the gas is 
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completely disrupted. If such solid cores survive tidal disruption of the extended gas 
envelope, they may very well represent the first rocky/icy planets or planetary embryos 
in a planetary system. 

Understanding the behavior of solids in gaseous disks is required for building a basic 
picture of planet formation, even in the context of a gravitationally unstable disk. As 
demonstrated here, fragments, although only modestly enriched, can be differentiated at 
birth. This leaves open the possibility of a hybrid scenario that we call core assist plus 
gas capture. In addition, it also suggests that clumps as massive as brown dwarfs could 
also start as differentiated objects, depending on the degree of the enhancement of solids 
within the spiral arms at the time of gas disk fragmentation. These formation channels 
should have observational consequences, which we only begin to touch upon here, but include 
metallicity and core size, as well as orbital locations if done in a statistical sense. Because 
the fragmentation process is strongly affected by dust characteristics, future studies should 
include grain growth models that self-consistently change dust opacities and capture, even 
if heuristically, the formation of solids with a distribution of sizes. Although challenging, 
such a model will greatly aid in determining the likelihood of fragmentation and what type 
of objects disk instability will typically form. 
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A. Drift Tests 

A difficulty witli our drag algoritlim is tliat it only includes the drag force when calculat- 
ing the new 5v. If the stopping time is much shorter than the hydrodynamic step, then the 
asymptotic behavior, e.g., radial drift velocities, could be poorly calculated. In this section, 
we present a simple test that addresses this issue. We find that the asymptotic solution 
is affected by our approximation, but for the regime studied in the simulations presented 
above, our approach is valid. 

Consider equation (1), but in the limit that the particles are subsonic, i.e., /e ~ 1. 
Take the gas velocity to be zero, so the velocity difference represents just the velocity of 
the particle. Assume that gravity, g, is constant and that the particle starts at rest to find 
Vd = gt'si^ — 6xp(— t/tg)), where t'^ = {n /SY^'^tg. The asymptotic drift velocity is simply 
gt'^. In the context of a disk, we can compute the asymptotic settling velocity of a particle 
by Vz = —Vl?zt'g^ assuming the stopping length is small. For the radial velocity, we follow 
Weidenschilling (1977) and consider the radial drift regime, where t'^ is small compared with 
the orbital time, as well as the perturbed orbit regime, where t'^ is long compared with the 
orbital time. For the radial drift regime, Vr = (^Ls ~ ^i-)^^s! ^^^ ^'^^ ^^e perturbed orbit 
regime, Vr = {^l^, - n]^)r/{t'^n]^). 

For this test, we envisage a disk that has T = 30 K and p = 10^^'^ g cm^^, everywhere. 
The central star is assumed to have a mass of 1.5 Mq, but the gas is fixed to orbit as if the 
star's mass were 1.35 Mq. The radial and vertical motions of the gas are assumed to be zero. 
For particles 10 cm in size or smaller, the particle is placed in the disk at an initial r = 100.2 
AU and z = 4.2 AU with an ^initial = ^i.smq- At sizes of a meter or larger, the particles 
need more time to adjust to the asymptotic solution, so they are placed at r = 200.2 AU 
and z = 4.2 AU. For the drift speeds, we measure Vr and Vz as the particles cross r = 95 
AU. We use a time step ~ 0.03 yr, which is similar to the hydrodynamic time step in SIMl 
and SIM2. 

We test particles that are 100 m, 10 m, 1 m, 10 cm, 1 cm, 1 mm, and 100 /im in radius, 
each with an internal particle density pa = 3 g cm~^. The results are shown in Figure 24. 
The general behavior is quite good through several orders of magnitude. The 100 /xm case 
does show stronger coupling than expected, and we attribute this to our assumption that 
only the drag force matters during the drag update. For the vertical settling speeds, we test 
particles sizes of 100/im, 1 mm, and 1 cm, for which we find the ratio between the actual and 
expected settling to be 0.94, 0.99, and 1.01, repsectively. Larger particles are not included 
in the settling test because they are not expected to reach the asymptotic solution. For 
example, the 10 cm particle has a t'^ ~ 200 yr, which is about 1/4 of the orbital period 
at r ~ 100 AU. For the disk parameters in these simulations, our drag algorithm seems to 
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be sufficiently accurate. The reader should also note that, in these regions of disks, 10 cm 
particles have the largest drift velocities, which is why we chose to explore this size scale. 
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SIM1, Initial Q 



SIM1, Initial Surface Density 




1000 




Fig. 1. — Initial Toomre Q and surface density S profiles for the ICs for SIMl. The particle 
surface density is multiplied by a factor of 100 in this plot. 
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Fig. 2. — Face-on surface density of the ICs for SIM2 . The colorbar shows log S(cgs). Spiral 
waves are clearly forming, and in about 1/2 of an orbit, the over density near x,y = 150, 150 
AU will form a dense spiral arm and fragment. 
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SIM2, Initial Q 



SiM2, initial Surface Density 





Fig. 3. — Initial Q and S profiles for the SIM2 ICs. The shallow infall profile causes mass 
to build up at large radii during disk formation, driving the outer disk toward extreme 
instability. The particle surface density is multiplied by a factor of 100 in this plot. 
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Fig. 4. — Toomre Q and S profiles just after injecting tlie disk onto tlie liigli pliysical 
resolution grid. The choppy Q profile is a result of the direct injection method, which 
creates a rough epicyclic frequency profile. Because this is only done once during the disk 
simulation, its effect is just to add extra noise at this particular snapshot. 
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Fig. 5. — Gas surface density snapshots for SIMlmm. Each square is 300 AU on a side, and 
the snapshots correspond to 340, 620, 780, 940, 1100, and 1260 yr going from left to right 
and top to bottom. The colorbar shows the logarithmic surface density in g cm~^. 
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Fig. 6. — Gas surface density snapshots for SIMlmu. Each square is 300 AU on a side, and 
the snapshots correspond to 620, 780, 940, 1100, 1260, 1410, 1570, 1720, and 1870 yr going 
from left to right and top to bottom. The colorbar shows the logarithmic surface density in 
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Fig. 7. — A plot of the gas temperature versus density values at all computational grid points 
for the 1 mm and the 1 /im opacity simulations corresponding to the 1260 yr snapshot. For 
the same temperature, the 1 /im case can reach higher densities, which eventually leads to 
fragmentation. The line called Cold Pressure, shows the thermodynamic temperature that 
would correspond to a given density for the given cold pressure. Clump formation occurs 
before this threshold is met, and cannot be responsible for altering fragmentation in the 1 
mm case. 
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Fig. 8. — Opacities from the 1mm and 1/im tables. For the temperatures in the spiral in 
these simulations, the 1 mm case is much more opaque. Only after clump formation do 
temperatures rise enough to make the 1mm opacities less opaque than the 1/im. 
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Fig. 9. — Gas surface density snapshots for SIMlmu (top) and SIMlmuNoRocks (bottom). 
Each square is 300 AU on a side, and the snapshots correspond to 1260, 1410, and 1570 yr 
going from left to right. The colorbar shows surface density in g cm~^. The evolutions are 
qualitatively similar between the disks, but differences, e.g., the number of fragments, are 
apparent. 
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l/.tm opacity with particles 




Fig. 10. — Gas surface density snapshots with rock positions superposed for SIMlmu. Each 
square is 300 AU on a side, and the snapshots correspond to 620, 780, 940, 1100, 1260, 1410, 
1570, 1720, and 1870 yr going from left to right and top to bottom. The colorbar shows the 
logarithmic surface density in g cm~^. 
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Fig. 11. — The logarithm of the rock-to-gas ratio in cyhndrical coordinates for several snap- 
shots of SIMlmu. The spiral arms show large surface density enrichment, while the interarm 
regions show significant solid depletion. The snapshots are for 620, 940, 1260, and 1870 
yr, from left to right and top to bottom. The left panel shows the ratio of the full surface 
densities, while the right panel shows the ratio only for the midplane cells. 
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Fig. 12. — Left: Azimuthally averaged surface density profiles of the rocks in SIMlmu. As 
spiral arms develop, the radial surface density profiles show strong radial concentrations of 
rocks and redistribution of material outward as well as inward. Right: The azimuthally 
averaged gas surface density shown with the rock surface density, multiplied by 100, at 1870 
yr in SIMlmu. Gas surface density peaks show corresponding rock density enhancements. 
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Fig. 13. — Gas surface density snapshots for SIM2. Each square is 500 AU on a side, and 
the snapshots correspond to 320, 450, 640, 800, 960, 1130, 1300, 1450, and 1620 yr going 
from left to right and top to bottom. The colorbar shows the logarithmic surface density in 
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SIM2 with particles 
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Fig. 14. — The same as in Figure 13, but with the particle distribution superimposed. 
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Fig. 15. — The logarithm of the rock-to-gas ratio in cyhndrical coordinates for the SIM2 
snapshots 640, 9600, 1300, and 1620 yr. The left panel shows the ratio of the full surface 
densities, while the right panel is only for the midplane cells. 
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Fig. 16. — Left: Azimuthally averaged surface density profiles of the rocks in SIM2. As 
in SIMlmu, spiral arms concentrate and redistribute the rocks. Right: The azimuthally 
averaged gas surface density shown along with the rock surface density, multiplied by 100, at 
1620 yr in SIM2. Gas surface density peaks show corresponding rock density enhancements. 
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Fig. 17. — Similar as Figure 14, but for SIM2kni. 
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Fig. 18. — The logarithm of the rock-to-gas ratio in cyhndrical coordinates for SIM2km. 
The snapshots are at similar times as shown in Figure 15. Only the ratio of the full surface 
densities is shown. The contours represent the surface density of the gas, demonstrating 
that the regions with enhanced solid-to-gas ratios are outside the gaseous spiral arms. In 
addition, although the clump itself is slightly depleted in solids, the regions immediately 
around the clump do have super solids-to-gas ratios. 
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Fig. 19. — Similar to Figure 15, but for a region zoomed in around the clump in SIM2. The 
left panel shows the logarithm of the rock-to-gas ratio for the entire surface density, while 
the right panel is just for the midplane cells. 



48 



SIM2:C1 



SIM2:C1 




R(AU) 



140 
120 
100 
80 
60 
40 
20 



1300 yr 
1420 yr 
1620 yr 




i^4»m 



10 



R(AU) 



12 



Fig. 20. — Cell values for the gas density and temperature as a function of distance from the 
density maximum in the SIM2:C1 clump. Several snapshots are shown. 
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Fig. 21. — The spiral arm in SIM2 just before fragmentation at 960 yr (left) and after 
fragmentation at the 1300 yr snapshot (right). Gas density is shown by the colorbar, while 
rock particles are given by black dots. The rocks that are considered to be part of the clump 
are shown with large, blue crosses. The same particles are shown in the 960 yr snapshot. 
Most of the mass comes from a 20-30 AU section of the spiral arm. 
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Fig. 22. — Solids associated with SIM2km:Cl at 1720 yr, and then traced back to an early 
stage of its orbit (light gray) and to later stages of the clump's orbit (black). Particles coming 
into the clump are not easily trapped, and can be scattered far from the planet's Hill sphere. 
The earliest snapshot shown is 1290 yr, and the latest is 2190 yr. 
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Fig. 23. — Cell values for the gas density and temperature as a function of distance from the 
density maximum in the SIMlmu:C2 clump. Several snapshots are shown. 
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Fig. 24. — Actual versus analytic asymptotic radial drift speeds for the solids. At very small 
particle sizes, the stopping times become comparable to the hydrodynamics time step, which 
causes more solid-gas coupling than expected. For the regime explored in the simulations 
presented in this paper, the algorithm is acceptable. 
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Table 1: Gas and rock mass for the clump that forms in SIM2, as well as one comparison with 
SIM2km. The disk radius r and azimuthal angle of the clump are given for each data set, 
and can be compared with Figure 13. Three different density thresholds are shown, where 
only mass above the threshold is considered. Full represents the entire mass of the clump, 
here for all gas that has p > 9 x 10~^^ g cm~^ and for T > 34 K. Half refers to all mass 
that is above half of the peak density in the clump, pmax, and > 90% is for only mass that is 
greater than 90% of the peak. This last threshold isolates the core conditions of the clump. 
The factor f^toG gives the ratio of rock mass to gas mass, and /cm = (/moG + 0.01)/0.02 
gives the total solids enrichment relative to the average nebula's value. For our conversion, 
we used IMj = 32OM0. We remind the reader that when we write "rocks," we are referring 
to a mixture of silicates and ices. 
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/moG 


/cnr 






Full 


8.1 




42 


0.016 


1.3 






Half 


3.5 




38 


0.034 


2.2 






>90% 


0.85 




32 


0.12 


6.4 






S1M2:C1 


Time = 1620 


yr 


r ~ 106 AU 


0-58° 


Pmax ~ 


= 1.7 X 10- 


-1 1 —9 

^^ g cm 


Density Threshold 


Gas (Mj) 




Rocks (Me) 


/moG 


/cnr 






Full 


8.5 




42 


0.015 


1.3 






Half 


3.5 




42 


0.037 


2.4 






>90% 


1.2 




38 


0.10 


5.5 






SlM2km:Cl 


Time = 1430 


yr 


r ~ 106 AU 


~ 346° 


Pmax ~ 


= 1.3 X 10" 


-11 g cm~^ 


Density Threshold 


Gas {Mj) 




Rocks (Me) 


/moG 


/cnr 






Full 


7.1 




7.2 


0.0032 


0.66 






Half 


3.7 




2.6 


0.0022 


0.61 






>90% 


1.3 




0.67 


0.0016 


0.58 
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Table 2: Same as Table 1, but for clumps in SIMlmu. Compare positions with Figure 6. 



SIMlmu:Cl 


Time = 1100 


yr 


r ~ 67 AU 


~ 330° 


Pmax 


= 2.3 X 


10- 


-11 


g cm 2 


Density Threshold 


Gas {Mj) 




Rocks (Me) 


/moG 


Jcnr 










Full 


8.2 




55 


0.021 


1.5 










Half 


3.1 




53 


0.052 


3.1 










>90% 


0.44 




38 


0.27 


14 










SIMlmu:C2 


Time = 1260 


yr 


r ~ 87 AU 


0-51° 


Pmax 


= 6.2 X 


10- 


-11 


g cm-2 


Density Threshold 


Gas (Mj) 




Rocks (Me) 


/moG 


Jcnr 










Full 


6.6 




38 


0.018 


1.4 










Half 


2.6 




27 


0.033 


2.1 










>90% 


0.66 




12 


0.059 


3.5 










SIMlmu:C2 


Time = 1410 


yr 


r ~ 91 AU 


(j) ~ 137° 


Pmax 


= 3.2 X 


10- 


-11 


g cm- 2 


Density Threshold 


Gas {Mj) 




Rocks (Me) 


/moG 


Jcnr 










Full 


11 




70 


0.020 


1.5 










Half 


3.4 




67 


0.062 


3.6 










>90% 


0.47 




48 


0.32 


16 










SIMlmu:C2 


Time = 1570 


yr 


r ~ 100 AU 


~ 214° 


Pmax 


= 4.6 X 


10- 


-11 


g cm"2 


Density Threshold 


Gas (Mj) 




Rocks (Me) 


/moG 


Jenr 










Full 


15 




83 


0.017 


1.4 










Half 


4.6 




77 


0.052 


3.1 










>90% 


0.57 




74 


0.40 


21 










SIMlmu:C2 


Time = 1720 


yr 


r ~ 104 AU 


~ 276° 


Pmax 


= 6.3x 


10- 


-11 


g cm- 2 


Density Threshold 


Gas {Mj) 




Rocks (Me) 


/moG 


Jcnr 










Full 


20 




110 


0.018 


1.4 










Half 


5.4 




86 


0.050 


3.0 










>90% 


0.53 




80 


0.47 


24 










SIMlmu:C2 


Time = 1870 


yr 


r ~ 100 AU 


~ 344° 


Pmax 


= 1.2 X 


10- 


-10 


gcm-2 


Density Threshold 


Gas (Mj) 




Rocks (Me) 


/moG 


Jcnr 










Full 


32 




170 


0.017 


1.3 










Half 


7.3 




170 


0.071 


4.1 










>90% 


1.4 




160 


0.35 


18 
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Table 3: £ 


lame as Table 1, but for clump C3 in SIMli 


mu. Compare posit 


ions with Figure 6. 


SIMlmu 


:C3 


Time = 1410 


yr 


r ~ 86 AU 


0~ 11° 


Pmax 


4.2 X 10-^1 g cm-2 


Density Threshold 


Gas {Mj) 




Rocks (Me) 


/moG 


Jcnr 




Full 




10 




54 


0.017 


1.4 




Half 




3.3 




48 


0.045 


2.8 




>90% 




0.29 




32 


0.34 


18 




SIMlmu 


:C3 


Time = 1570 


yr 


r ~ 82 AU 


0-97° 


Pmax 


4.7 X 10"" g cm-2 


Density Threshold 


Gas {Mj) 




Rocks (Me) 


/moG 


Jcnr 




Full 




11 




61 


0.017 


1.4 




Half 




3.4 




58 


0.053 


3.1 




>90% 




0.61 




51 


0.26 


14 




SIMlmu 


:C3 


Time = 1720 


yr 


r ~ 77 AU 


~ 198° 


Pmax 


3.0 X 10"" g cm-2 


Density Threshold 


Gas {Mj) 




Rocks (Me) 


/moG 


Jcnr 




Full 




10 




64 


0.020 


1.5 




Half 




3.2 




61 


0.059 


3.5 




>90% 




0.36 




61 


0.53 


27 




SIMlmu 


:C3 


Time = 1870 


yr 


r ~ 86 AU 


~ 297° 


Pmax 


4.3 X 10-" g cm-2 


Density Threshold 


Gas {Mj) 




Rocks (Me) 


/moG 


Jcnr 




Full 




11 




64 


0.018 


1.4 




Half 




3.7 




64 


0.054 


3.2 




>90% 




0.64 




66 


0.33 


17 





